##########################################################################################

library(data.table)
library(optparse)
library(ggplot2)
library(Seurat)
library(ggsci)
library(ArchR)
library(ggplotify)
library(magrittr)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--atac_file"), type = "character"),
    make_option(c("--scriptPath"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## atac文件
    atac_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## 既往研究整理的代码
    scriptPath <- "~/20231121_singleMuti/scripts/scScalpChromatin"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
atac_file <- opt$atac_file
scriptPath <- opt$scriptPath
out_path <- opt$out_path

dir.create( paste0(out_path , "/markerGenes") , recursive = T)

###########################################################################################

a <- load(rna_file)
DefaultAssay(scrnat) <- "RNA"
## scrnat

b <- load(atac_file)
## testis_combined_peak_combineRNA

##########################################################################################
## 已发表文献写好的脚本
source(paste0(scriptPath, "/plotting_config.R"))
source(paste0(scriptPath, "/misc_helpers.R"))
source(paste0(scriptPath, "/matrix_helpers.R"))
source(paste0(scriptPath, "/archr_helpers.R"))
source(paste0(scriptPath, "/GO_wrappers.R"))

###########################################################################################

featureSets <- list(
    "SSC" = c("GFRA1","UTF1"), 
    "Differenting&Differented SPG" = c("KIT","DMRT1","STRA8"), 
    "Leptotene" = c("SPO11","SYCP3"), 
    "Zygotene" = c("SYCP3"), 
    "Patchytene" = c("OVOL1","OVOL2"),
    "Diplotene" = c("NME8"), 
    "Early stage of spermatids" = c("C9orf116","ACR"),
    "Round&ElongateS.tids" = c("TXNDC8","TXNDC2"), 
    "Sperm" = c("TNP1","PRM2"), 
    "Leydig cells" = c("DLK1"), 
    "Sertoli cells" = c("AMH","SOX9"), 
    "Myoid cells" = c("MYH11","ACTA2"), 
    "Pericytes" = c("NOTCH3"),
    "Macrophages" = c("CD68","CD14"),
    "Endothelial cells" = c("VWF","PECAM1"), 
    "NKT cells" = c("NKG7","FGFBP2")
)

grp_order = c("SSC", "Differenting&Differented SPG", "Leptotene",
"Zygotene", "Patchytene", "Diplotene",
"Early stage of spermatids", "Round&ElongateS.tids", "Sperm",
"Leydig cells", "Myoid cells", "Pericytes",
"Sertoli cells", "Endothelial cells", "NKT cells", "Macrophages")

gene_order <- c("GFRA1", "UTF1", "KIT",
"DMRT1", "STRA8", "SPO11",
"SYCP3", "OVOL1", "OVOL2",
"NME8", "C9orf116", "ACR",
"TXNDC8", "TXNDC2", "TNP1",
"PRM2", "DLK1", "MYH11",
"ACTA2", "NOTCH3", "AMH",
"SOX9", "VWF", "PECAM1",
"NKG7", "FGFBP2", "CD68",
"CD14")

gene_order <- gene_order[length(gene_order):1]

###########################################################################################
## RNA气泡图
namedClustAspect <- 1.6

count_mat <- GetAssayData(object=scrnat,assay="RNA",layer="counts")
count_mat <- as.matrix(count_mat)
avgPctMat <- avgAndPctExpressed(count_mat, scrnat$cell_type, feature_normalize=TRUE, min_pct=0)

subGenes <- featureSets %>% do.call("c",.)
avgPctMat <- avgPctMat[avgPctMat$feature %in% subGenes,]
avgPctMat <- avgPctMat[avgPctMat$grp %in% grp_order,]

avgPctMat$pctExpr[avgPctMat$pctExpr < 5] <- 0
        
p5 <- dotPlot(
    avgPctMat, xcol="grp", ycol="feature", color_col="avgExpr", size_col="pctExpr", xorder=grp_order, yorder=gene_order, cmap=cmaps_BOR$sunrise, aspectRatio=namedClustAspect) + scale_size_continuous(range=c(1,4)) + 
    theme(axis.text.x=element_text(size=8), axis.text.y = element_text(size=8,face="italic"), panel.border = element_rect(size = 0.2), axis.ticks = element_blank(), legend.title = element_text(size = 8), legend.text = element_text(size = 8), legend.key.size = unit(0.3, "cm")) +
    xlab(NULL) + ylab(NULL) + guides(
    fill = guide_legend(title=""),
    colour = guide_colorbar(title="Relative expression", override.aes = list(size = 3, linetype = 2, shape = 16)),
    size = guide_legend(title="Percentage expressed"))

out_file <- paste0( out_path , "/RNA_NamedClust_markers_dot_plot_scalp.pdf" )
ggsave(out_file , p5 , width=5,height=6)

###########################################################################################
## ATAC气泡图
GSM_se <- getMatrixFromProject(testis_combined_peak_combineRNA, useMatrix="GeneScoreMatrix")
GSM_mat <- assays(GSM_se)$GeneScoreMatrix
rownames(GSM_mat) <- rowData(GSM_se)$name

## 转化为数值矩阵
GSM_mat_num <- apply(GSM_mat , 1 , as.numeric)
rownames(GSM_mat_num) <- colnames(GSM_mat)
GSM_mat_num <- t(GSM_mat_num)

avgPctMat <- avgAndPctExpressed(GSM_mat_num[,getCellNames(testis_combined_peak_combineRNA)], testis_combined_peak_combineRNA$cell_type, feature_normalize=TRUE, min_pct=0)

subGenes <- featureSets %>% do.call("c",.)
avgPctMat <- avgPctMat[avgPctMat$feature %in% subGenes,]
avgPctMat <- avgPctMat[avgPctMat$grp %in% grp_order,]

avgPctMat$pctExpr[avgPctMat$pctExpr < 5] <- 0

p6 <- dotPlot(
    avgPctMat, xcol="grp", ycol="feature", color_col="avgExpr", size_col="pctExpr", xorder=grp_order, yorder=gene_order, cmap=cmaps_BOR$horizonExtra, aspectRatio=namedClustAspect) + scale_size_continuous(range=c(1,4)) + 
    theme(axis.text.x=element_text(size=8), axis.text.y = element_text(size=8,face="italic"), panel.border = element_rect(size = 0.2), axis.ticks = element_blank(), legend.title = element_text(size = 8), legend.text = element_text(size = 8), legend.key.size = unit(0.3, "cm")) +
    xlab(NULL) + ylab(NULL) + guides(
    fill = guide_legend(title=""),
    colour = guide_colorbar(title="Relative gene activity", override.aes = list(size = 3, linetype = 2, shape = 16)),
    size = guide_legend(title="Percentage expressed"))

out_file <- paste0( out_path , "/GSM_NamedClust_markers_dot_plot_scalp.pdf" )
ggsave(out_file , p6 , width=5,height=6)